function detect_opening(I, center, radius)
    cx = center(1);
    cy = center(2);
    range = cy-radius:cy+floor(radius/2);
    
    r = I(:,:,2);
    gfilter = fspecial('gaussian',10,1);
    r = imfilter(r,gfilter);
    
    leftpart = r(:,5:cx-radius-10);
    rightpart = r(:,cx+radius+10:end-5);
    
    [~,gl] = gradient(double(leftpart));
    [~,gr] = gradient(double(rightpart));
    
    points_left = zeros(size(leftpart));
    pl = zeros(range,2);
    pr = zeros(range,2);
    points_right = zeros(size(rightpart));
    
    figure(1);
    subplot(2,1,1); hold on; imshow(leftpart);
    subplot(2,1,2); hold on; imshow(rightpart);
    for i=range
        % Get the local minima
        currentvalue = leftpart(i,:);
        avg = sum(currentvalue)/length(currentvalue);
        currentvalue = smooth(double(currentvalue),30);
        [~,~,ymin,imin2] = extrema(currentvalue);
        
        imin2(ymin > avg) = -1;
        
        imin2 = sort(imin2,'descend');
        %points_left(i,imin2(1)) = 1;
        pl(i,1) = imin2(2);
        
        subplot(2,1,1);
        plot(imin2(2),i,'rx');
        
        currentvalue = rightpart(i,:);
        currentvalue = smooth(double(currentvalue),15);
        [~,~,~,imin2] = extrema(currentvalue);
        
        % Note: it doesn't have descend here
        imin2 = sort(imin2);
        %points_right(i,imin2(2)) = 1;
        pr(i,1) = imin2(2);
        subplot(2,1,2);
        plot(imin2(2),i,'rx');
        
        % Get the brightness jump
        subplot(2,1,1);
        currentvalue = gl(i,:);
        currentvalue = smooth(double(currentvalue),20);
        [~,imax2] = max(currentvalue);
        points_left(i,imax2) = 1;
        pl(i,2) = imax2;
        plot(imax2,i,'bx');
        
        subplot(2,1,2);
        currentvalue = gr(i,:);
        currentvalue = smooth(double(currentvalue),15);
        [~,imax2] = max(currentvalue);
        points_right(i,imax2) = 1;
        pr(i,2) = imax2;
        plot(imax2,i,'bx');
    end
    
    figure(1);
    [H,T,R] = hough(points_left,'RhoResolution',5,'Theta',0:89.5);
    P  = houghpeaks(H,1,'threshold',ceil(0.3*max(H(:))));
    lines = houghlines(points_left,T,R,P,'FillGap',5,'MinLength',7);
    xy = [lines(1).point1; lines(1).point2];
    yleft = interp1(xy(:,2),xy(:,1),range,'linear','extrap');
    subplot(2,1,1); plot(yleft,range);
    
    [H,T,R] = hough(points_right,'RhoResolution',5,'Theta',-89.5:0);
    P  = houghpeaks(H,1,'threshold',ceil(0.3*max(H(:))));
    lines = houghlines(points_right,T,R,P,'FillGap',5,'MinLength',7);
    %subplot(2,1,2);
    xy = [lines(1).point1; lines(1).point2];
    yright = interp1(xy(:,2),xy(:,1),range,'linear','extrap');
    subplot(2,1,2); plot(yright,range);
    
    yl = [];
    xl = [];
    a = abs(pl(range,2)-yleft');
    figure(5); imshow(r); hold on;
    for i=range
        if a(i-(cy-radius)+1) < 10
            yl = [yl i];
            xl = [xl pl(i,2)];
            plot(pl(i,2),i,'rx');
        end
    end
    
    figure(5);
    yr = [];
    xr = [];
    a = abs(pr(range,2)-yright');
    for i=range
        if a(i-(cy-radius)+1) < 4
            yr = [yr i];
            xr = [xr pr(i,2)+cx+radius];
            plot(pr(i,2)+cx+radius,i,'rx');
        end
    end
    
    x = [xl xr];
    y = [yl yr];
    
    f = polyfit(x,y,3);
    curve = polyval(f,1:size(r,2));
    plot(curve);
end